### Code reads CESOP survey and examine questions about knowledge of programs/policies
### Produces the figure and  table in Appendix H 
### CESOP's usage policy does not allow us to distribute the microdata
### Data is publicly available in CESOP's website: https://www.cesop.unicamp.br
### Look for survey 04683
### Please contact authors at natalia.bueno@emory.edu or cesar.zucco@fgv.br 
### if you are unsuccessful in accessing the data via CESOP

### Saves tables in .tex format to subfolder "Tables"
### Saves figures in .pdf format to subfolder "Figures"

rm(list=ls())
library(questionr)
library(car)
library(foreign)
library(xtable)
library(MASS)
### Some functions ######

nsnr <- function(x){  
  if(class(x)=="numeric"|class(x)=="integer"){
    ifelse(is.element(x,c(88,98,99)),NA,x)
  }else{
    x <- gsub("^\\s?-|–\\s","",x)
    ifelse(is.element(tolower(as.character(x)),
                      c("ns/nr","ns","nr","nsa","não sabe/ não respondeu","não respondeu","não sabe","não respondeu/ recusa", "não sabe (esp.)", "não respondeu (esp.)")),NA,as.character(x))}
}

my.regs <- function(x,xs,a=.95){
  require(plotrix)
  baseline <- mean(x$model[,1])
  points <- summary(x)$coef[,1]
  ciwd <- summary(x)$coef[,2]*qnorm(a)
  output <- list(baseline=baseline,est=cbind(points,ciwd))
  return(output)
}

my.plot <- function(xx,the.xs=1,sym.col=T){
  #workhouse function for plot using the same yscale for all quesitons
  require(plotrix)
  parties <- nrow(xx$est)
  xs <- the.xs+seq(-.275,.275,len=parties)  
  if(sym.col==T){
    plotCI(xs,xx$est[,1],uiw=xx$est[,2]
           ,col=c(1,gray(0.25),gray(0.25),gray(0.25),1),lwd=1.5,sfrac=0.005
           ,pch=NA,ylab="",xlab="",yaxt="n",xaxt="n",bty="n",add=T)
    points(xs,xx$est[,1],pch=21:(21+parties)
           ,bg=c(1,rep("white",parties-2),1)
           ,col=c(1,rep(gray(0.45),parties-2),1))}
  if(sym.col==F){
    plotCI(xs,xx$est[,1]+xx$baseline,uiw=xx$est[,2]
           ,col=c(1,gray(0.25),gray(0.25),gray(0.25),gray(0.25)),lwd=1.5,sfrac=0.005
           ,pch=NA,ylab="",xlab="",yaxt="n",xaxt="n",bty="n",add=T)
    points(xs,xx$est[,1]+xx$baseline,pch=21:(21+parties)
           ,bg=c(1,rep("white",parties-1))
           ,col=c(1,rep(gray(0.25),parties-1)))}		
  
  segments(x0=xs[1],x1=xs[parties],y0=xx$baseline,y1=xx$baseline,lty=2,lwd=1.5)
}

my.baseplot<-function(x,ymin=0.5,ymax=1,my.axis=T){#background empty plot
  plot(c(.5,x+.5),c(ymin,ymax),type="n",xlab="",ylab="",xaxt="n",yaxt="n",bty="n")
  for(i in 1:x){polygon(x=c(i-.5,i+.5,i+.5,i-.5),
                        y=c(-1,-1,11,11),border=NA,
                        col=ifelse(i/2==round(i/2),gray(0.85),"white"))}
  if(my.axis){axis(side=2,at=seq(ymin,ymax,by=1),labels=seq(ymin,ymax,by=1))}
}

the.plot  <- function(x,ymin=3,ymax=9,my.axis=T){#creates the plot with same scale, calls other ploting fcts
  my.baseplot(length(x),ymin=ymin,ymax=ymax,my.axis=my.axis)
  xx<- lapply(x,my.regs)
  for(i in 1:length(x)){my.plot(xx[[i]],i)}
}

### Read and recode some background variables 

d <-  read.spss('DataPrivate/04683.SAV'
                ,to.data.frame=T
                ,add.undeclared.levels = c("sort")
                ,duplicated.value.labels = c("append"))
names(d)<-tolower(names(d))

### Recode income and a few others to predict missing income

d$fam.income <- factor(car::recode(d$rend2,"'ATÉ 1'='00.01';
                                    'MAIS DE 1 A 2'='01.02';
                                   'MAIS DE 2 A 5'='02.05';
                                   'MAIS DE 5 A 10'='05.10';
                                   c('MAIS DE 10 A 20','MAIS DE 20')='10.'; 
                                   'NÃO TEM RENDIMENTO PESSOAL'='00.01';
                                   'NÃO RESPONDEU'=NA"),ordered=T)
d$ind.income <- factor(car::recode(d$rend1,"'ATÉ 1'='00.01';
                                   'MAIS DE 1 A 2'='01.02';
                                   'MAIS DE 2 A 5'='02.05';
                                   'MAIS DE 5 A 10'='05.10';
                                   c('MAIS DE 10 A 20','MAIS DE 20')='10.'; 
                                   'NÃO TEM RENDIMENTO PESSOAL'='00.01';
                                   'NÃO RESPONDEU'=NA"),ordered=T)


d$pent <- d$p45=="Assembléia de Deus"|
          d$p45=="Universal do Reino de Deus"|
          d$p45=="Deus é Amor"|
          d$p45=="Evangelho Quadrangular"|
          d$p45=="Igreja Internacional da Graça"|
          d$p45=="Renascer em Cristo"|
          d$p45=="Sara nossa terra"|
          d$p45=="Outras Evangélicas específicas"|
          d$p45=="Evangélica - Não sabe especificar"

d$white <- d$p44=="Branca"|d$p44=="Amarela"

### Do individuals or family members participate in BF?
### Individuals could list up to 8 social programs (eight variables)
### And for family members they could list up to 5
### So, we merge the 8 into one, and the 5 into one indicating participatin in BF only

d$bf.family <- F
bf.a <- c(grep("Bolsa Família",d$p46a01),
            grep("Bolsa Família",d$p46a02),
            grep("Bolsa Família",d$p46a03),
            grep("Bolsa Família",d$p46a04),
            grep("Bolsa Família",d$p46a05),
            grep("Bolsa Família",d$p46a06),
            grep("Bolsa Família",d$p46a07),
            grep("Bolsa Família",d$p46a08))
bf.b <- c(grep("Bolsa Família",d$p46b01),
            grep("Bolsa Família",d$p46b02),
            grep("Bolsa Família",d$p46b03),
            grep("Bolsa Família",d$p46b04),
            grep("Bolsa Família",d$p46b05))

d$bf.family [union(bf.a,bf.b)] <- T
d$superior <- d$inst=="SUP. INC."|d$inst=="SUP. COMP."
d$analfa <- d$inst=="ANALF."
d$N.NE <- d$reg=="NORTE"|d$reg=="NORDESTE"
d$capital <- d$condx=="CAPITAL"

# Predicting regression (not using individual income because there's NAs there too)
pred.reg <- polr(fam.income~pent+white+bf.family+N.NE+capital+analfa+superior, d, method = c("logistic"))
summary(pred.reg)

# Predict missing income observations
d$fam.income.imp <- d$fam.income
d$fam.income.imp[is.na(d$fam.income)]<-   predict(pred.reg,newdata = d[is.na(d$fam.income),])
d$fam.income.imp <- factor(d$fam.income.imp,ordered=F) #remover ordering

### Do individuals or family members participate in MCMV?
### Individuals could list up to 8 social programs (eight variables)
### And for family members they could list up to 5
### So, we merge the 8 into one, and the 5 into one indicating participatin in MCMV only

d$mcmv.individual <- d$mcmv.family <- F
mcmv.a <- c(grep("Minha Casa",d$p46a01),
            grep("Minha Casa",d$p46a02),
            grep("Minha Casa",d$p46a03),
            grep("Minha Casa",d$p46a04),
            grep("Minha Casa",d$p46a05),
            grep("Minha Casa",d$p46a06),
            grep("Minha Casa",d$p46a07),
            grep("Minha Casa",d$p46a08))
mcmv.b <- c(grep("Minha Casa",d$p46b01),
            grep("Minha Casa",d$p46b02),
            grep("Minha Casa",d$p46b03),
            grep("Minha Casa",d$p46b04),
            grep("Minha Casa",d$p46b05))

d$mcmv.individual[mcmv.a] <- T
d$mcmv.family[union(mcmv.a,mcmv.b)] <- T


# Pergunta espontanea
# O(a) sr(a) se lembra de alguma ação, programa ou obra do Governo Federal que esteja em andamento neste momento? (CASO SIM) De qual ou quais ações, programas ou obras o(a) sr(a) se lembra? Mais alguma? 
# Qualmais? (ESPONTÂNEA – ATÉ 3 MENÇÕES - CASO NÃO ENCONTRE NA LISTA OU NÃO CONSIGA ENCAIXAR A RESPOSTA, ANOTE DA FORMA MAIS COMPLETA POSSÍVEL – CASO DEIXE DE CITAR ALGUMA MENÇÃO, MARQUE OPÇÃO ‘NÃO RESPONDEU’) 
# Até três respostas por pessoa

d$prog1 <- car::recode(d$p1101,"'Bolsa Família'='bf';
                                'Minha Casa Minha Vida'='mcmv';
                                'Casas populares/ Projetos voltados à habitação'='mcmv';
                                'Mais Médicos'='mm';
                                'ENEM - Exame Nacional do Ensino Médio'='enem';
                                'ProUni - Programa Universidade para Todos '='prouni';
                                'Fies - Fundo do Financiamento Estudantil'='fies';
                                'Pronatec - Programa Nacional de Acesso ao Ensino Técnico'='pronatec';
                                'Farmácia Popular'='farmpop';
                                c('Não lembra','Não sabe','Não respondeu')='nsnr';
                                else='other'")
d$prog2 <- car::recode(d$p1102,"'Bolsa Família'='bf';
                                'Minha Casa Minha Vida'='mcmv';
                       'Casas populares/ Projetos voltados à habitação'='mcmv';
                       'Mais Médicos'='mm';
                       'ENEM - Exame Nacional do Ensino Médio'='enem';
                       'ProUni - Programa Universidade para Todos '='prouni';
                       'Fies - Fundo do Financiamento Estudantil'='fies';
                       'Pronatec - Programa Nacional de Acesso ao Ensino Técnico'='pronatec';
                       'Farmácia Popular'='farmpop';
                      c('Não lembra','Não sabe','Não respondeu',NA)='nsnr';
                       else='other'")
d$prog3 <- car::recode(d$p1103,"'Bolsa Família'='bf';
                                'Minha Casa Minha Vida'='mcmv';
                       'Casas populares/ Projetos voltados à habitação'='mcmv';
                       'Mais Médicos'='mm';
                       'ENEM - Exame Nacional do Ensino Médio'='enem';
                       'ProUni - Programa Universidade para Todos '='prouni';
                       'Fies - Fundo do Financiamento Estudantil'='fies';
                       'Pronatec - Programa Nacional de Acesso ao Ensino Técnico'='pronatec';
                       'Farmácia Popular'='farmpop';
                       c('Não lembra','Não sabe','Não respondeu',NA)='nsnr';
                       else='other'")

cited <- function(program,data=d){
  d$prog1==program|d$prog2==program|d$prog3==program
}

all.cited <- data.frame(
          mcmv=cited("mcmv"),
          bf=cited("bf"),
          prouni = cited("prouni"),
          pronatec = cited("pronatec"),
          enem = cited("enem"),
          farmpop = cited("farmpop"),
          fies = cited("fies"),
          mm = cited("mm"),
          other = cited("other"),
          none = d$prog1=="nsnr"&d$prog2=="nsnr"&d$prog3=="nsnr")

most.cited.espontanea <- apply(all.cited,2,sum)
most.cited.espontanea <- data.frame(N=most.cited.espontanea,                                 
                          pct=round(apply(all.cited,2,sum)/nrow(d)*100,2))


# Pergunta estimulada (P1201-P1227,P1297-P1299)
# O(a) sr(a) tomou conhecimento, mesmo que de ouvir falar de alguma destas ações,
# programas ou obras do Governo Federal que estejam em andamento neste momento? 
# (CASO SIM) Qual ou quais ações,

est.all.cited <-apply(d[,grep("^p12",names(d))],2,function(x){sum(x=="Sim")})
sort(est.all.cited,decreasing=T) #the top seven are the same we found in the "spontaneous"
d<- rename.variable(d,"p1201","bf") 
d<- rename.variable(d,"p1202","pronatec")
d<- rename.variable(d,"p1203","mcmv")
d<- rename.variable(d,"p1208","mm")
d<- rename.variable(d,"p1209","prouni")
d<- rename.variable(d,"p1210","fies")
d<- rename.variable(d,"p1214","farmpop")
d<- rename.variable(d,"p1225","enem")

know.vars <- c("mcmv","bf","prouni","pronatec"
               ,"enem","farmpop","fies","mm")
d[,know.vars] <- sapply(d[,know.vars], function(x){x=="Sim"})

d$other <- d$p1204=="Sim"|d$p1205=="Sim"|d$p1206=="Sim"|
  d$p1207=="Sim"|d$p1211=="Sim"|d$p1212=="Sim"|d$p1213=="Sim"|
  d$p1215=="Sim"|d$p1216=="Sim"|d$p1217=="Sim"|d$p1218=="Sim"|
  d$p1219=="Sim"|d$p1220=="Sim"|d$p1221=="Sim"|d$p1222=="Sim"|
  d$p1223=="Sim"|d$p1224=="Sim"|d$p1226=="Sim"|d$p1227=="Sim"

d$none <- d$p1297=="Sim"|d$p1298=="Sim"|d$p1299=="Sim"

selected <- d[,c("mcmv","bf","prouni","pronatec","enem","farmpop","fies","mm","other","none")]
most.cited.estimulada <- apply(selected,2,sum)
most.cited.estimulada <- data.frame(N=most.cited.estimulada ,                                 
                                    pct=round(most.cited.estimulada/nrow(d)*100,2))
                                    
cited.table <- merge(most.cited.espontanea,most.cited.estimulada
                        ,by=0,all=T,suffixes=c(" Esp"," Est"))


## Evaluation of social programs
## Programs are grouped into three separate questions P13, P14, P15
## Here we select only the eight main programs, listed in the previous quesitons
d<- rename.variable(d,"p1501","ev.bf") 
d<- rename.variable(d,"p1502","ev.mcmv")
d<- rename.variable(d,"p1402","ev.farmpop")
d<- rename.variable(d,"p1301","ev.enem")
d<- rename.variable(d,"p1304","ev.pronatec")
d<- rename.variable(d,"p1401","ev.mm")
d<- rename.variable(d,"p1303","ev.fies")
d<- rename.variable(d,"p1302","ev.prouni")
eval.vars <- c("ev.mcmv","ev.bf","ev.prouni","ev.pronatec"
               ,"ev.enem","ev.farmpop","ev.fies","ev.mm")
d[,eval.vars] <- sapply(d[,eval.vars], function(x){as.numeric(nsnr(x))})

# Average evaluation
selected <- d[,eval.vars]
ev.mean <- apply(selected,2,mean,na.rm=T)

# Dispersion
ev.sd <- apply(selected,2,sd,na.rm=T)

# Add eval info to table with knowledge
ev.table <- merge(ev.mean,ev.sd,by=0)
ev.table$Row.names <- gsub("ev.","",ev.table$Row.names)
the.table <- merge(cited.table,ev.table,by="Row.names",all.x=T)

# Compare evaluation of MCMV and BF
t.test(na.omit(d$ev.mcmv),na.omit(d$ev.bf),alternative="two.sided")

# Compara evaluation of MCMV and BF among non-beneficiaries of the respective program
mean(d$ev.mcmv[which(d$mcmv.family==F)],na.rm=T)
mean(d$ev.bf[which(d$bf.family==F)],na.rm=T)
t.test(na.omit(d$ev.mcmv[which(d$mcmv.family==F)]),na.omit(d$ev.bf[which(d$bf.family==F)]),alternative="two.sided")

# Evaluation by income level(plot all)
all.lm <- function(i){
  the.formula <- formula(paste(eval.vars[i],"~fam.income.imp-1"))
  out <- lm(the.formula,data=d)
}
all.regs <- lapply(1:length(eval.vars),all.lm)
names(all.regs)<- gsub("ev.","",eval.vars)

pdf(file="Figures/fig-H1.pdf",width=12,height=6)
par(mar=c(4,4,2,1))
the.plot(all.regs[order(ev.mean,decreasing = T)],ymin=4,ymax=8)
mtext(side=2,"Evaluation of Program",line=2.5)
legend(x="bottom",
       legend=c("0 to 1 MW","1 to 2 MW","2 to 5 MW","5 to 10 MW","> 10 MW","Average"),
       pch=c(21:25,NA),col=1,pt.bg=c(1,"white","white","white",1,NA),
       bty="n",lty=c(rep(1,5),2),horiz=F,ncol=2
       ,cex=1,inset=-0.175,xpd=NA
) 
axis(side=3,at=c(1:8),tick=F,line=-3
     ,labels=c("MCMV","Bolsa\nFamilia","PROUNI","PRONATEC","ENEM","Farmacia\nPopular","FIES","Mais\nMédicos")[order(ev.mean,decreasing = T)]
     ,cex.axis=1)
dev.off()

## Evaluation by income level 
# (store difference between highest and lowest bracket)
all.lm2 <- function(i){
  the.formula <- formula(paste(eval.vars[i],"~fam.income.imp"))
  out <- lm(the.formula,data=d)
}
all.regs2 <- lapply(1:length(eval.vars),all.lm2)
all.regs2 <- do.call("rbind",lapply(all.regs2,function(x){summary(x)$coef[5,]}))
rownames(all.regs2)<- gsub("ev.","",eval.vars)

#add info to table for paper
table.difinc <- data.frame(Row.names=rownames(all.regs2),
                    DifInc=all.regs2[,"Estimate"])
the.table <- merge(the.table,table.difinc,by="Row.names",all.x=T)
x.the.table <- xtable(the.table[order(the.table$'pct Esp',decreasing = T),c(1,3,5,6,7,8)],digits=1)
caption(x.the.table) <- "Knowledge and Evaluation of Social Programs (2015)"
label(x.the.table) <- "tab:progs"
print(x.the.table,include.rownames=F
      ,caption.placement = "top"
      ,hline.after = c(0,nrow(x.the.table))
      ,file = "Tables/tab-h1.tex")
